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For classical billiards we suggest that a matrix of action or length of trajectories in conjunction with statistical 
measures, level spacing distribution and spectral rigidity, can be used to distinguish chaotic from integrable 
systems. As examples of 2D chaotic billiards we consider the Bunimovich stadium billiard and the Sinai billiard. 
In the level spacing distribution and spectral rigidity we find GOE behavior consistent with predictions from 
random matrix theory. The length matrix elements follow a distribution close to a Gaussian. We present for the 
Lorentz gas model a mathematical result stating that the length of trajectories starting from random points on 
the billiard boundary converge in distribution to a universal Gaussian in the limit when the number of collisions 
tends to infinity. As example of a 2D integrable billiard we consider the rectangular billiard. We find very rigid 
behavior with strongly correlated spectra similar to a Dirac comb. These findings present numerical evidence for 
universality in level spacing fluctuations to hold in classically integrable systems and in classically fully chaotic 
systems. Finally, we address the question if universality holds in chaotic potential models. 

PACS numbers: 05.40.-a, 05.45.-a 



I. INTRODUCTION 



The idea to model apparently disordered spectra, like 
those of of heavy nuclei, using random matrices was 
suggested around 1950 by Mehta and Wigner. They 
showed that random matrices of Gaussian orthogo- 
nal ensembles (GOE) generate a Wigner-type nearest- 
neighbor level spacing (NNS) distribution |1]. In a 
classic paper Bohigas, Giannoni and Schmit (BGS) for- 
mulated a conjecture |2], stating that time-reversal in- 
variant quantum systems with classically fully chaotic 
(ergodic) counterpart have universality properties given 
by random matrix theory (RMT). Experiments in nu- 
clear physics, for example, have shown that spectra 
originating from different heavy nuclei give the same 
Wignerian energy level spacing distribution |3]. Uni- 
versality properties in quantum chaos of bound sys- 
tems, i.e. quantum systems with a fully chaotic clas- 
sical counter part, have now been demonstrated in many 
experiments, computational models and in theoretical 
studies idllllHlTlli. Theoretical support of the BGS 
conjecture came from the semiclassical theory of spec- 
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tral rigidity by Berry |H, [13], who showed that univer- 
sal behavior in the energy level statistics is due to long 
classical orbits. Sieber and Richter 1 1 1] investigated the 
role of correlations between classical orbits. The semi- 
classical theory has been further developed by Miiller et 
al. 1 12, 13, 14]. In ref. (l^ they presented the 'core of 
a proof of the BGS conjecture, which, by proving ar- 
guments previously used by Berry (91], show that in the 
semi-classical limit the periodic classical orbits deter- 
mine the universal fluctuations of quantum energy lev- 
els. Further refinements have been made by Keating and 
Miiller ill. 

The study of classical strongly chaotic systems 
(Anosov systems) has revealed that central limit theo- 
rems (CTL) hold iHIiTiN, 19, 20, 21]. This has been 
proven for the 2D periodic Lorentz gas with finite hori- 
zon. The first step of proof was given by Bunimovich 
and Sinai (l^ and was completed by Bunimovich, Sinai 
and Chernov ifTsI] . At macroscopic times, such deter- 
ministic chaotic system converges to Brownian motion, 
i.e. behaves like a random system |[i3,[i8llSl2il]- This 
is also supported by the existence of an average diffu- 
sion coefficient |[T6[[T8[|22|] . showing the diffusive char- 
acter of such chaotic system. 

Do classical fully chaotic systems also exhibit uni- 
versality properties? This question was addressed by 
Argaman et al. ll23h . who showed that there is univer- 
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sal behavior in 2-point correlation functions of actions 
of classical periodic orbits of classically chaotic sys- 
tems. As examples they considered the deformed cat- 
map and the baker-map. This has been elaborated fur- 
ther in a number of studies by Dittes et al. [i24i] . Au- 
rich and Sieber iH, Cohen et al. lEl Tanner jllll, 
Sano ll28n. Primack and Smilansky 123], Sieber and 
Richter lflTIl and Smilansky and Verdene ||3QIi . Arga- 
man et al. started from the assumption that spectral 
fluctuations of chaotic quantum systems follow the pre- 
dictions of RMT and they derived a universal expres- 
sion for classical correlation functions of periodic or- 
bits via Gutzwiller's semi-classical trace formula. They 
concluded "The real challenge, though, is to find out 
whether these action correlations can be explained on a 
completely classical level". 

This was answered recently by Laprise et al. ll3lll 
who found universal behavior in classical 2D billiards 
by looking at fluctuations in spectra of classical ac- 
tion/length matrices from bouncing trajectories. In par- 
ticular, they considered the Lima^on/Robnik family of 
billiards, which interpolates between the chaotic car- 
diod billiard and the integrable circle billiard. For the 
cardioid billiard, the level spacing distribution P{s) and 
spectral rigidity A3 (L) were found to be consistent with 
the prediction of GOE behavior. For the interpolating 
case close to the circle, the behavior approached a Pois- 
sonian distribution. The circle billiard itself was found 
to be very rigid and strongly correlated and yielded 
P{s) oc 5(5' — 1). The jump in behavior at the transition 
to the circle is associated with the corresponding change 
in the symmetry group. 

This article extends the results of Ref. ll3l|] in the 
following directions: (i) We consider the 2D rectan- 
gular billiard as another example of an integrable bil- 
liard. Compared to the circular billiard, this billiard has 
lower symmetry (no group property). Nevertheless, it 
displays strong spectral correlation and rigidity like the 
circular billiard, (ii) We present numerical studies for 
other chaotic 2D billiards, like the Sinai-billiard and 
the Bunimovich stadium billiard, (iii) In order to un- 
derstand the observed universal behavior in chaotic bil- 
liards, we present arguments linking such behavior to 
CTLs, diffusive and random walk behavior, which hold, 
e.g., for the periodic Lorentz gas model with finite hori- 
zon. In particular, we present a mathematically rigorous 
result on the distribution of length of trajectories, (iv) 
We address the question if the observed universal be- 
havior holds more generally beyond billiards, e.g., for 
chaotic systems governed by a potential. We discuss 
some tentative results for optical billiards with variable 



refraction coefficients. 

The answers we found can be summarized as follows: 
For the 2D Sinai billiard and 2D Bunimovich stadium 
billiard the level spacing distribution P{s) and spectral 
rigidity A3 (L) are consistent with predictions of RMT 
(GOE behavior), i.e. show universal behavior. This 
behavior is statistically the same as the one observed 
in quantum chaos (obtained from energy level spacing 
distributions). The implication of these findings is that 
RMT not only represents well the statistical fluctuation 
properties of the energy spectrum of chaotic quantum 
systems, but also those of the length spectrum of chaotic 
classical systems. Moreover, statistical fluctuations ob- 
tained from spectra of action/length matrices clearly dis- 
tinguish chaotic from integrable systems. 



II. LENGTH AND ACTION OF TRAJECTORIES 

In classical systems, chaos information is encoded 
in trajectories. According to the Alekseev-Brudno the- 
orem ll32ll the temporal length t is related via the 
Kolmogorov- Sinai K entropy to the information I{t) in 
the segment of trajectory, 

limt^ooI{t)/t = K . (1) 

This motivates us to look at the length of trajectories and 
its fluctuation properties. Let L{x{t)^x{t)^t) denote the 
Langrangien of a system, let x^^^^ (^) denote a solution 
of the Euler-Lagrange equations of motion. Let 

Z = ^[^->]= rdtL{x{t),x{t),t)\,^,rraJ , 

A[/'-j]= rds\,^,.aj (2) 

denote the action and length, respectively, for such 
trajectory. We choose a finite set of discrete points 
{xi,X2, . . . ,Xiv}. For any pair of boundary points Xi^xj 
we compute a classical trajectory Xx^^lj, connecting 
those points. We suggest the construction of an action 
matrix and a length matrix, 

E/j = S[Xj^.^^^] , 

A,7=A[<%.]. (3) 

Both matrices are viable for statistical analysis of clas- 
sical chaos. In the case of billiard systems, we consider 
trajectories where the billiard particle moves with con- 
stant velocity u and constant kinetic energy E. Then the 
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action and the length matrix are essentially equivalent, 
^ij = ^^ij' (4) 



III. INTEGRABLE BILLIARD 

If one considers integrable quantum systems and an- 
alyzes them in terms of NNS energy level spacing dis- 
tribution and spectral rigidity, then in most cases one 
finds a Poissonian distribution P{s) which is reflected 
also in the behavior of A3(L). However, there are ex- 
amples of integrable quantum systems in 2D where the 
levels are correlated and the level spacing distribution is 
not given by a Poissonian. Such cases are called non- 
generic. Berry and Tabor |[33ll noted as example uncou- 
pled quantum oscillators in 2D. Casati et al. ll34ll and 
later Seligman and Verbaarschot |[35|] showed that also 
the integrable 2D quantum rectangular incommensurate 
billiard does not give an uncorrected Poissonian level 
spacing distribution. In a recent study of 2D quantum 
harmonic oscillators Chakrabarti and Hu found in 
the case of uncoupled oscillators a level spacing distri- 
bution displaying a 5(5' — 1 ) peak plus some background. 
They measured the spectral correlation via the spec- 
tral rigidity A3 (L) and observed a flat curve for L > 4 
saturating at A3(L)|^^^ = 0.17. Similar results were 
found in the case with weak harmonic coupling, yield- 
ing A3(L)|^^^ = 0.12 and correlation C = 0.997. This is 
close to the rigidity of a picket fence (Dirac comb) of 
equally spaced levels with A3(L) = 0.083. They con- 
clude that the 2D quantum oscillator system is highly 
correlated at short and long-range, is regular and very 
rigid. 

Laprise et al. ll3lll considered the classical integrable 
circle billiard and constructed a length matrix from 
classical trajectories between boundary points located 
evenly on the billiard wall. In the spectral rigidity they 
found saturation for large L at A3(L) = 0.11 which is 
close to the rigidity of a Dirac comb (A3(L) = 0.083). In 
conjunction with a correlation coefficient of C = 0.999 
and a level spacing distribution P{s) oc — 1) (Dirac 
comb) this indicates high correlation at short and long 
range and very rigid behavior. 

A. Rectangular Billiard 

As example of an integrable classical billiard we con- 
sider the 2D rectangular billiard, shown in Fig. |[T1. The 
shape is determined by the parameters a,b, which were 
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FIG. 1 : 2D Rectangular billiard. Horizontal length a, vertical 
length b. Trajectories go from boundary points Xi to xj located 
in the upper right quarter of the billiard wall. 




FIG. 2: Rectangular billiard. Top panel: Number of trajec- 
tories versus number of rebounds. Bottom panel: Error vs. 
number of rebounds. 

chosen to be a = ^/5 and b = I. The boundary points 
{xi, . . . ,Xiv}, located in the upper right quarter of the 
billiard wall, are distributed regularly, with perimeter 
spacing given by 

d{xn^uxn) = 2^js^^i^ ' 1,...,A^-1 . (5) 

For a given pair of boundary points, we find the behav- 
ior of the number of trajectories versus the number of 
rebounds, shown in Fig. O (top). The error behavior of 
trajectories as function of the number of rebounds has 
been obtained by taking into account Ntraj = 20 trajec- 
tories corresponding to different initial conditions. The 
results are presented in Fig. O (bottom), showing that 
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FIG. 3: Rectangular billiard. Top panel: Distribution of 
length matrix elements P(A). Middle panel: Spectral rigid- 
ity A3 (L) of length spectrum. Bottom panel: Level spacing 
distribution P{s) of length spectrum. 



the error is stable in the regime < Nreb < 100- This is 
in contrast to an exponential increase found in chaotic 
billiards. 

The symmetry of the rectangular billiard is mirror 
symmetry under reflection of the x-axes and of the y- 
axes (with origin at center of rectangle). The sym- 
metry shows up in the shape of trajectories. For ex- 
ample, a trajectory (1) going from starting point xn to 
endpoint xn-i bouncing once from the lower bound- 
ary wall has the same shape as the trajectory (2) go- 
ing from starting point to endpoint xa^_2 bouncing 
once from the lower boundary wall. Trajectory (2) is 
obtained by a translation of trajectory (1) in x-direction 



by the amount of perimeter spacing (Eq. O. However, 
such discrete translations do not form a group (they 
would form a group if one would consider the rectan- 
gular billiard with periodic boundary conditions at all 
billiard boundary walls, which would map the billiard 
on a torus, and the symmetry group would then be a 
rotation group). The length matrix is not exactly a cir- 
cular matrix either, however, shares with circular ma- 
trices the property of having a number of different ma- 
trix elements equal to or less than the rank of the ma- 
trix. This property implies strong correlations between 
length matrix elements, which translates to strong cor- 
relations of eigenvalues of the length matrix. To sum 
up, although there is no exact symmetry group, there 
are residues of a "broken symmetry group of discrete 
translations", which imply strong correlations among 
length matrix elements and among eigenvalues of the 
length matrix. We expect that this will manifest itself in 
the statistical behavior in the level spacing distribution 
P{s) and the spectral rigidity As{L). The results corre- 
sponding to the parameters N = 11 and Nreb = 55 are 
shown in Fig. O. Fig. O (top) represents the distri- 
bution P{A) of length matrix elements. Its shape looks 
quite different from the near-gaussian shape for chaotic 
billiards (see below). Fig. (Si (middle) shows the spec- 
tral rigidity A3 (L) of the spectrum of the length matrix. 
It rapidly saturates to a value A3(L)|^^^ = 0.088, which 
is close to the value A3{L)\Dimc = 1/12 ^ 0.083 of an 
ideal Dirac comb. We have also computed the correla- 
tion coefficient to obtain C = 0.996. This is consistent 
also with the NNS level distribution P{s) of the length 
matrix shown in Fig. O (bottom) where one observes 
P{s) oc 8(5' — 1) indicating the behavior of a Dirac comb. 
The rectangular billiard turns out to be highly correlated 
at short and long range, to be very rigid and to be regu- 
lar. 

In the classically integrable 2D rectangular billiard 
we found levels to be strongly correlated leading to a 
Dirac comb (picket fence) distribution. This is possibly 
evidence for universal behavior in the integrable case. 
Comparing the behavior of the rectangular billiard with 
the circular billiard ll3lll , we observe that they differ in 
their symmetry properties. In the circular billiard hop- 
ping from one boundary point to its neighbor stands for 
a group operation. The corresponding operation in the 
rectangular billiard has no group property. However, 
the resulting strong correlation and spectral rigidity are 
found to be very similar for both billiards. 

Compared to the behavior of integrable quantum sys- 
tems, which is typically of Poissonian nature with cor- 
relation coefficient C = 0, the behavior of classically in- 
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tegrable system is completely opposite, with correlation 
coefficient C ^ I. 



IV. STOCHASTIC BILLIARD SYSTEM 



A. Random Walk 
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We look at the statistical behavior of A for many trajec- 
tories, having in common the same number Nreb of col- 
lisions. For statistical purpose it is meaningful to hold 
Nreb fixed, which means that we keep fixed the average 
length of trajectory 

{A) = {Nreb^l){Ax) , (8) 

and A fluctuates around (A) . When we plot a histogram 
of length P{A), we obtain a Gaussian distribution, as 
shown in Fig.lH. This result holds provided that we 
consider macroscopic times, i.e. a large number of col- 
lisions. Such behavior can be understood from the CTL: 
In between collisions, any length segment Ak = Axk of 
the straight trajecory is a random number drawn from 
some distribton, say P{x). The total length A then is 
a sum of random numbers. Under rather general condi- 
tions on the distribution P{x), the CTL says that the sum 
of random variables becomes a Gaussian in the limit of 
large Nreb- The physical meaning is this: The length 



FIG. 4: Random walk for Nreb = 1000. Histogram of length 
matrix elements. Gaussian behavior is predicted by the CTL. 

As pointed out in sect. Jl fully chaotic deterministic 
systems (like the Lorentz gas) may display behavior of 
a random system, like Brownian motion, diffusive char- 
acter and validity of CTLs. Thus, before discussing 
chaotic billiards, we first take a look at a stochastic 
billiard. The collision and diffusion of molecules in 
gaseous or liquid matter can be described by Brow- 
nian motion, which can be modeled by an open ran- 
dom billiard. In the limit of free path length Ax ^ 0, 
At ^0 and {Ax)^/Ax = const, the model describes dif- 
fusion and is mathematically described by the following 
differential equation, having a probabilistic interpreta- 
tion irlllsllli, 



a,p=Z)(V)2p 



(6) 



Let us consider the simple model of a massive particle, 
representing, e.g., an atom or a molecule, which from 
time to time collides with other atoms /molecules. Af- 
ter each collision, the particle changes its direction and 
velocity. We assume that the new direction and velocity 
are randomly distributed. In between collisions, the par- 
ticle moves freely, maintaining direction and velocity. 
The length of trajectory of the random walker having 
executed Nreb collisions is given by 



A = ^k= X Ay^ . 

k=l k=l 
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FIG. 5: Gaussian random matrix of dimension 5000. Top 
panel: Distribution of Gaussian random matrix elements. 
Middle panel: Spectral level spacing distribution P{s). Bot- 
tom panel: Spectral rigidity A^{L). 

of path A of a random walker between any given start- 
ing point Xin and any end point Xfi follows a Gaussian 
distribution, if we keep fixed the number Nreb of colli- 
sions, and if Nreb is sufficiently large. If we construct 
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FIG. 6: Uniform random matrix of dimension 5000. Top 
panel: Distribution of uniform random matrix elements. Mid- 
dle panel: Spectral level spacing distribution P{s). Bottom 
panel: Spectral rigidity A3(L). 



a length matrix A^^, where initial and end points are 
taken from a given set of points {xi, . . . ,Xiv}, then all 
matrix elements follow a Gaussian distribution (same 
variance and mean). Because Afj is real and symmet- 
ric, A represents a matrix from a Gaussian orthogonal 
ensemble (GOE). Random matrix theory predicts that 
the level spacing distribution of the eigenvalues from a 
GOE matrix follows a Wigner distribution. 

This is similar to the statistical approach to quantum 
chaos, where the level spacing distributions of energy 
eigenvalues from a quantum Hamiltonian are predicted 
to behave according to random matrix theory (i.e. like 
GOE in the case of a time-reversal symmetric system). 
According to the BGS-conjecture, this holds when the 
classical counter part of the quantum system is fully 
chaotic. 

Although the random walker is a purely random sys- 
tem, and not a deterministic chaotic system, such ob- 
served behavior of the classical length of the random 
walker in analogy to the behavior of Hamiltonian matrix 
elements in chaotic quantum systems led us to postu- 
late the following hypothesis: With respect to statistical 



fluctuations, the length matrix may play the same role 
in classical chaos as the Hamilton matrix in quantum 
chaos. In the following we investigate this hypothesis. 

Mathematical note. The BGS-conjecture does not 
state that the matrix elements itself of a quantum Hamil- 
tonian are distributed like a GOE ensemble. The con- 
jecture rather says only that the statistical fluctuations 
of the eigenvalue spacings obtained from the quantum 
Hamiltonian are the same as those from a GOE ensem- 
ble, giving a Wignerian distribution. In other words, it is 
possible that the matrix elements of the quantum Hamil- 
tonian are distributed quite differently from a Gaussian, 
but nevertheless its level spacing distribution is Wigne- 
rian. A mathematical example is given in Figs. II5I6I . 
Fig. O shows the behavior of Gaussian random matrix, 
yielding a level spacing distribution P{s) of Wignerian 
type (GOE) and also a spectral rigidity A3 (L) consistent 
with GOE. However, within statistical errors the same 
behavior in P{s) and A3(L) is obtained from a uniform 
random matrix (non-GOE) as shown in Fig. (61. 

Such a situation, where the distribution of matrix el- 
ements is not GOE, but the level fluctuation statics is 
GOE, occurs in nuclear physics. An example is the 
distribution Hamiltonian matrix elements obtained from 
nuclear shell model calculations 1 6] . In this model, there 
are vanishing Hamiltonian matrix elements. This im- 
plies that the number of independent matrix elements 
is much smaller than in a random matrix of the same 
size. However, Brody et al. ||4Q|i could show that the 2- 
body residual interaction in the shell model yield matrix 
elements of random character following a Gaussian dis- 
tribution. In particular, they showed that spectral fluc- 
tuation properties from such ensembles with orthogonal 
symmetry are identical to those from GOE. That implies 
that GOE is meaningful to predict spectral fluctuation 
properties of nuclei governed by 2-body interactions. 



V. CHAOTIC BILLIARDS 

A. Sinai Billiard 

For general closed 2D billiards, the mean free path 
length A free in between two collisions is given by the 
billiard geometry via ll4lll 



^free — 



m 



(9) 



where \Q\ stands for the bilUard area and \dQ\ for the 
length of the wall. Let us now consider a particle of 
mass m moving in the 2D Sinai billiard (see Fig. Q). 
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FIG. 7: Geometry of 2D Sinai billiard. Parameters: R = I, 
d = 5. Trajectories go from boundary points Xi to xj located 
in a segment of one eighth of the wall of the billiard. 
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FIG. 9: 2D Sinai billiard. Nreb = 15. Top panel: Distribution 
of length matrix elements. Middle panel: Level spacing distri- 
bution from length matrix A. Bottom panel: Spectral rigidity. 
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FIG. 8: 2D Sinai billiard. Top panel: Number of trajectories 
for a given number of rebounds. Bottom panel: Relative error 
as function of number of rebounds. 



The classical Sinai billiard system is known to be fully 
chaotic 1I42I1 . The billiard is symmetric under mirror op- 
eration about X- and y-axes, and about the diagonals of 
angles 7r/4 and 37r/4, respectively, passing through the 
center. In order not to mix different symmetry classes, 
we consider the billiard with boundary points xi , . . . ,Xiv 
located in one eighth of the exterior billiard wall (see 
Fig- El)- It is important to make sure that the chaotic 
behavior comes from the dynamics of the system and 
not from random location of boundary points. Thus we 
have located the boundary points in a regular way, with 



N-\ 



, ^= 1,...,A^-1 . (10) 



The rule of dynamics is free motion in the interior region 
and perfectly elastic specular collision at the central disc 
and the exterior square wall. According to Eqs. l8l9l the 
length of trajectories fluctuates around 



(A) = (A^,,^ + 1) 



(11) 



We have classified trajectories according to the num- 
ber Nreb of rebounds. For a given number of rebounds 
Nreb and a given pair of boundary points x^, xj, we de- 
termined the trajectories by starting from Xi with an an- 
gle a = —71/2 (measured with respect to the normal to 
the boundary at point Xi). We gradually increased a 
and searched for a value, say ai, such that the trajec- 
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tory arrived at the endpoint xj. We continued to in- 
crease a until we found a value a2, such that another 
trajectory arrived at xj. In such a way we obtained a 
set of Ntraj trajectories, corresponding to starting angles 
av,V = 1, . . . ^ Ntraj (the value of Ntraj depends on the 
particular pair of boundary points). We repeated this 
for all combination of boundary points and thus created 
a set of length matrices A^^^ v = l,2,3, — A global 
characteristic feature of chaos is encoded in the num- 
ber of classical trajectories. For the Sinai billiard we 
found that the number A/^^; of trajectories averaged over 
all pairs of boundary points increases as function of the 
number of rebounds Nreb like an exponential, 

Ntraj = exp[a A^,^^] - 1 , a = 0.81 ±0.03 , (12) 

shown in Fig. (HI (top panel). It turns out that such 
exponential behavior in chaotic billiards is clearly dis- 
tinct from the behavior in integrable billiards, where 
the number of trajectories increases linearly with the 
number of rebounds (see rectangular billiard). By hav- 
ing chosen a value of Nreb and keeping it fixed, we 
have generated an ensemble of length matrices { A^ | v = 
1, . . . , Ntraj} corresponding to all possible trajectory in- 
dices V. In general one expects that chaotic behavior 
develops with increasing Nreb- 

In order to make sure that the chaotic behavior is not 
due to numericl noise, we estimated the numerical error 
\xfi —^in\ by following a trajectory from a given start- 
ing point Xin until it carried out Nreb rebounds, then fol- 
lowing the trajectory in time reversed direction to ar- 
rive after another Nreb at some x/i. We observe that the 
error behavior has two regimes (Fig. O bottom). For 
Nreb < 48 it shows an exponential increase, following 
on average the rule 

£ = exp[PA^,^^-41] , P = 0.76 ±0.03 forA^,^^ <48 . 

(13) 

For Nreb > 48 saturation is reached (the error is in the 
order of one). Such exponential behavior is related to 
the largest positive Lyapunov exponent of the system 
and thus represents a global characteristic of a chaotic 
system. Like the exponential behavior of number of tra- 
jectories (Eq. [12]), also the exponential error behavior 
distinguishes between chaotic billiards and integrable 
bilHards (Fig. d). 

From Fig. (81 (bottom) we find a relative error of 
about 10~^^ for Nreb < 15. After unfolding the spec- 
tra we retain 8 significant digits. We used the technique 
of gaussian broadening (vl to smooth the raw spec- 
trum and unfold the spectrum. In this method there is 
free parameter, which we tuned to reproduce the auto- 



correlation coefficient (C = — 0.27 for Wigner distribu- 
tion d, 113]). For a given trajectory index v we com- 
puted a level spacing distribution. Afterwards we super- 
imposed the spectra corresponding to different trajec- 
tory indices V. For Nreb = 15 the number of trajectories 
is quite large (Ntraj ~ 10^). In order to limit compu- 
tational cost, we have restricted ourselves by consider- 
ing smaller subsets of length matrices, in the follow- 
ing way. Above we have discussed for a given number 
of rebounds Nreb and a given pair of boundary points 
Xi, Xj how we determined the starting angles (mea- 
sured from the normal to the boundary line at the bound- 
ary point) which corresponds to a trajectory v. That 
search was carried out in the range — 7l/2<a<7l/2. 
Such subsets have been constructed by choosing (i) a re- 
stricted range of starting angles, e.g. —7z/3 < a < 7c/2, 
or — 7l/4 < a < 7l/2, etc. and (ii) also introducing an 
upper bound ntraj Ntraj) on the number of trajecto- 
ries and retaining only the trajectories corresponding to 
the starting angles ai , . . . , oc^^^^y . We repeated this for all 
combination of boundary points and thus created differ- 
ent subsets of length matrices A^^^ , with v = 1 , . . . , ntraj- 

It turns out that ntraj = 20 gives a reasonable good 
statistics. We verified that the resulting distributions 
are independent of the particular subset of trajectories 
within statistical errors. The resulting level spacing 
distribution P{s) of length eigenvalues and the spec- 
tral rigidity A3(L), corresponding to A/^ = 50 boundary 
points, Nreb = 15 rebounds and ntraj = 20 trajectories, 
are shown in Fig. (9l (middle and bottom panel). The 
results show a Wigner distribution for P{s), and are con- 
sistent with GOE behavior also for A^i^L). As in the 
quantum Sinai billiard, the level spacing distribution has 
been shown to give a Wigner distribution [2] . 

The histogram of the length matrix elements itself is 
shown in Fig. (51 (top panel). The distribution looks 
close to a Gaussian. Is it a pure Gaussian? This ques- 
tion is physically relevant for the following reason: If 
the distribution P{A) is a Gaussian, then the matrix ele- 
ments obey GOE statistics. Then RMT [1] implies that 
the level fluctuation statistics P{s) and the spectral rigid- 
ity A3 (L) follow GOE statistics (which seems consisent 
with the numerical results shown in Fig. (9l (middle and 
bottom panel). We pointed out above that the random 
walk model gives a Gaussian distribution for the his- 
togram of the length matrix elements P{A) . The random 
walk model is mathematically a Markov chain, i.e., the 
length of each piece of straight path in between subse- 
quent collisions is given by a random number. Two sub- 
sequent random numbers are statistically independent. 
On the other hand, the chaotic Sinai billiard is a deter- 
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ministic system, and the length of two subsequent pieces 
of straight trajectory are not independent. Hence it is not 
plausible that the distribution /^(A) of length matrix el- 
ements for the Sinai billiard is given by a Gaussian. In 
that case the question is: Why does the level flucuation 
statistics follow GOE? We will discuss this at the end of 
this section. 



B. Bunimovich Stadium Billiard 




FIG. 10: 2D stadium billiard. Trajectories go from bound- 
ary points Xi to Xj located in a segment of one quarter of the 
billiard wall. 



[Linear regression 
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Let us consider the 2D Bunimovich stadium billiard 
with semi-axes a and b (Fig. ifTOl ). The billiard is known 
to be fully chaotic [45] . The billiard is symmetric under 
mirror operation about x- and y-axes. In order not to 
mix different symmetry classes, we consider the billiard 
with boundary points xi , . . . ,Xiv located in one quarter of 
the billiard wall. The quarter perimeter has the length 
b -\-7l{b — a)/2. The nodes Xn are equi-distantly dis- 
tributed on the quarter perimeter of the stadium, with 
perimeter spacing given by 

d{xn^uxn) = [b^n{b-a)/2]/{N-l) n=l,...,N-l . 

(14) 

The rule of dynamics is free motion in the interior re- 
gion and perfectly elastic specular collision at the exte- 
rior square wall. For a given pair of boundary points, 
the number of trajectories Ntraj connecting these points 
increases on average with the number of rebounds ex- 
ponentially, according to 

Ntraj = exp[a (A^,,^ - 1)] , a = 1.01 ±0.09 , (15) 

shown in Fig. ifTTIl (top). Such behavior is qualitatively 
similar to that found in the Sinai billiard. We also mea- 
sured the numerical error following the method used in 
the Sinai billiard. Likewise, we found a regime of expo- 
nential behavior followed by a regime of saturation de- 
picted in Fig. ifm (bottom). On average the exponential 



FIG. 11: 2D Stadium billiard. Stadium parameters a/b = 2.2. 
Statistical average over 15 trajectories. Top panel: Number of 
trajectories Nsoi as function of number of rebounds Nyeb- Bot- 
tom panel: Relative error 8 as function of number of rebounds 

Nreb- 



increase is given by 

8 = exp[PA^,^^-33.1] , P= 1.13 ±0.05 forA^,^^<28. 

(16) 

We have classified trajectories by the the number of 
rebounds Nreb^ which was kept fixed. In this way, 
we generate an ensemble of length matrices {ly |v = 
1, . . . J Ntraj} corresponding to all possible trajectory in- 
dices V. For each trajectory index v we computed a 
level spacing distribution. Afterwards we superimposed 
the spectra corresponding coming from a randomly cho- 
sen considerably smaller subset of trajectory indices 
V = 1 , . . . , ritraj (like we did for the Sinai billiard, see 
above). The resulting level spacing distribution P{s) of 
length eigenvalues and the spectral rigidity A3(L) are 
shown in Fig. {T2\ (middle and bottom panel). The re- 
sults are consistent with a Wigner distribution, i.e. GOE 
behavior. The histogram of the length matrix elements 
P{A) shown in Fig. |[T2ll (top) looks close to a Gaussian. 

In order to see if the observed Wigner distribution in 
the level spacing distribution depends on the statistical 
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FIG. 12: 2D stadium billiard. Average over trajectories. Ge- 
ometry parameters a/b = 2.2. Boundary points on quarter of 
wall N = 40. Top panel: Distribution of length matrix ele- 
ments. Nreb = 30. Middle panel: Level spacing distribution 
P{s) from length matrix A. Nreb = 12. Average over ritraj = 18 
trajectories. Bottom panel: Spectral rigidity A3 (L) (same pa- 
rameters as middle panel). 



method of averaging over several trajectories (i.e. length 
matrices A^) we have tested an alternative, where we av- 
eraged over different stadium geometries. All stadium 
billiards are fully chaotic, and if the Wigner behavior 
is universal, this should show up for any stadium ge- 
ometry. Superimposing spectra corresponding to differ- 
ent stadium shapes should improve statistics. In their 
paper on quantum chaos in nuclear physics, Bohigas 



et al. ||46|] have analyzed spectra of a series of differ- 
ent heavy nuclei (but having the same quantum num- 
bers) and showed that their level fluctuation statistics 
agrees with the prediction of RMT. These different nu- 
clei correspond to different potentials. In analogy to 
that we considered here length spectra corresponding 
to different stadium shapes. The stadium wall repre- 
sents a curve where the potential jumps from zero to 
infinity. Hence different stadium shapes correspond to 
different potentials. By superimposing spectra from dif- 
ferent stadium shapes we obtained results displayed in 
Fig. |[T3l . Within statistical errors the results are equiv- 
alent to those obtained by superimposing spectra from 
different trajectories (Fig. |[T2ll ). 



C. Universality in chaotic billiards 

The numerical experiments with chaotic billiards in- 
vestigated above led us to the following observations: 

(i) For a given pair of boundary points and a given num- 
ber of rebounds Nreb, the number of classical trajectories 
Ntraj shows exponential behavior as function of Nreb- 

(ii) For a given pair of boundary points, the numeri- 
cal error 8 shows exponential behavior as function of 
the number of rebounds Nreb- (iii) The NNS level spac- 
ing fluctuation distribution P{s) obtained from eigenval- 
ues of the length matrix A shows a Wignerian distribu- 
tion. This is consistent with GOE behavior predicted 
by RMT. (iv) The spectral rigidity A3 (L) obtained from 
the spectrum of length matrix A is also found consistent 
with GOE behavior predicted by RMT. (v) The distribu- 
tion P{A) of length matrix elements Aij is found to be 
close to a Gaussian. 

The leading Gaussian behavior in the distribtion P{A) 
is possibly universal. How can we understand such be- 
havior of P{A) in chaotic billiards? In the case of the 
random walk model the distribution P{A) was shown to 
be a Gaussian, being a direct consequence of the CTL. 
Let us now consider the chaotic Sinai billiard, which 
is equivalent to a system where the billiard ball moves 
in an open system of equal spherical discs located on a 
2D rectangular regular grid, known as the 2D (diluted) 
Lorentz gas model. The billiard ball alternates between 
free motion and collisions with the circular discs. There 
are two types of 2D periodic Lorentz gases (or Sinai bil- 
liards on a torus): One has a finite horizon, where free 
paths between collisions are bounded (the scatterers are 
so dense that they block every direction of motion. An 
example on a triangular lattice is given in Ref. I22i1 ). 
The other type has an infinite horizon, where the parti- 
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FIG. 13: 2D stadium billiard. Average over shapes of sta- 
dium. Top panel: Distribution of length matrix elements. 
Nreb = 12. Middle panel: Level spacing distribution from 
length matrix A. Bottom panel: Spectral rigidity A3. 



cle can move indefinitely without collision. The Sinai 
billiard investigated above (see Fig.||7l) belongs to this 
class. The CTL, existence of finite diffusion coefficient 
and convergence to Brownian motion was proven only 
for the periodic Lorentz gas with finite horizon |[l6[ [T8ll. 
In the case of infinite horizon, the CTL has not been 
proven. Moreover, in this case the diffusion coefficient 
is infinite and there is no convergence to Brownian mo- 
tion iil. 

Actually, for the periodic Lorentz gas with finite hori- 
zon it can be shown rigorously that the distribution of 



length of trajectories becomes a Gaussian distribution 
in the limit of many bounces. This holds when the ini- 
tial points of trajectories are distributed randomly on the 
billiard boundary. Then Chernov and Markarian ll48ll 
prove the following result 




(17) 

for all — 00 < ^ < 00. Here x(x) denotes the return func- 
tion and 

tn = x(x) + x{Fx) + • • • + x(F"-^x) (18) 

is the time of the n-th collision, is the variance of the 
random variable depending on the observable /, which 
here is the return function x(x). F : M ^ M is the col- 
lision map and M is collision space (phase space of bil- 
liard map). Thus the distribution of the travel time until 
the n-th collision obeys the CTL, and hence converges 
in distribution to a normal, i.e. Gaussian distribution in 
the limit where the number of collisions goes to infin- 
ity. In the case where the billiard particle moves with 
constant speed u, this means also the length of trajec- 
tory In = utn obeys the CTL, and the length distribution 
P{ln) tends to a normal Gaussian distribution for n^oo. 
This result seems to support our numerical findings 
of a (near) Gaussian distribution of length of trajecto- 
ries. However, the scenario where the above mathe- 
matical result holds and the scenario of our numerical 
study differ in two respects, namely in the distribution 
of boundary points and in the horizon of billiard. For the 
purpose of statistical analysis in terms of RMT we are 
interested in the distribution P(A^), which corresponds 
to the case where the billiard particle starts from and 
arrives at boundary positions taken from a discrete set 
of boundary points, which is different from an initial 
random set of boundary points. Nevertheless, it seems 
plausible that a similar mathematical result may hold 
also for the discrete set of initial and final boundary 
points, such that P(A^) tends to a Gaussian. Second, 
the mathematical result holds for the case of finite hori- 
zon, while the Sinai billiard in our numerical study has 
an infinite horizon. At a first glance, it looks puzzling in 
absence of strong mathematical results (for infinite hori- 
zon) that our numerical results nevertheless show (near) 
Gaussian behavior for the length distribution. A possi- 
ble explanation is this: The choice of the discrete set of 
regularly distributed boundary points selects preferen- 
tially trajectories which often bounce at the central disc, 
while trajectories with infinite horizon (no bounces at 
disc) are avoided. This is supported also by numerical 
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evidence that the distribution P(A^) of the Sinai billiard 
(Fig.O top) and that of the rectangular biUiard viewed 
as Sinai bilHard with central disc of radius zero (Fig.O 
top) are quite different. 

For the stadium billiard, Balint and Gouezel [l49i] have 
shown that the CTL also holds. They proved that the 
limit distribution (n oo) is Gaussian, but the scaling 
factor is not given by the standard expression y^, but 
rather by ^nlogn. As consequence, however, there 
is no finite diffusion coefficient and no convergence to 
Brownian motion. Balint and Gouezel prove that the 
distribution of length of trajectories, tends to a normal 
distribution in the limit when the number of collisions 
n goes to infinity. Like in the Lorentz gas, this corre- 
sponds to a random distribution of initial points on the 
boundary. This is consistent with our numerical results 
of the stadium billiard which show in P{A) small devi- 
ations from a Gaussian comparable in magnitude with 
those in the Sinai model. They are likely due to the 
small number of bounces. 

Let us suppose that P{An) does tend to a Gaussian. 
As a consequence, in this limit (and after suitable nor- 
malisation) the length matrix A itself will become a 
GOE matrix 1 5] . This is sufficient to guarantee in this 
limit that the level fluctuations P{s) and spectral rigid- 
ity A3(L) obtained from the matrix A show GOE be- 
havior 111]. This does establish universal behavior in 
the limit of many bounces in the case of the Lorentz 
gas/Sinai billiard as well as the Bunimovich stadium bil- 
liard. 



VI. IS THERE UNIVERSALITY IN CLASSICALLY 
CHAOTIC POTENTIAL SYSTEMS? 

In quantum chaos, universal behavior has been ob- 
served experimentally in billiard systems llsoll as well 
as in potential systems, like the hydrogen atom in a 
strong magnetic field 1 5 1 ] . Thus one may expect uni- 
versal behavior to show up also in classically chaotic 
potential systems. In order to investigate this, in a first 
step, we have considered an optical billiard. This is 
given by a stadium shaped billiard. Several smaller sta- 
dium billiards are located inside (similar to russian pup- 
pets). The resulting geometry looks like tracks in a track 
and field stadium. To each "track" we assign an opti- 
cal medium of a given constant refraction coefficient, 
in a way such that the refraction coefficient monoton- 
ically increases in steps going from the inner field to 
the outer track. For the outer stadium boundary we im- 
pose total reflection. A light ray moving in such sce- 



nario will undergo a change of direction when cross- 
ing a "track" boundary, following Snell-Descartes' law. 
Likewise it will change its speed, also following Snell- 
Descartes 'law. Such change in direction and speed is 
similar to the effect of the presence of a potential. In 
a numerical experiment analogous to that of the Buni- 
movich stadium billiard we made a statistical analysis 
of trajectories. Again we found universal behavior con- 
sistent with the prediction of GOE. More details will be 
given elsewhere. We consider this finding as a tentative 
result hinting towards the possibility of universal behav- 
ior in classically chaotic potential systems. 



VII. SUMMARY 

This paper is about classical chaos occuring widely 
in nature, for example in astrophysics, meteorology and 
dynamics of the atmosphere, fluid and ocean dynam- 
ics, climate change, chemical reactions, biology, phys- 
iology, neuroscience, or medicine. We have suggested 
to extend random matrix theory, used in chaotic quan- 
tum systems, to classically chaotic and integrable sys- 
tems. We have studied fully chaotic as well as integrable 
billiards and used a statistical description based on the 
length of trajectories to discriminate chaotic versus in- 
tegrable behavior. 
Results: 

(i) In chaotic billiards (stadium and Sinai billiard) the 
NNS distribution P{s) as well as A3(L) obtained from 
the length matrix show GOE behavior, i.e., they are 
universal. The implication is that RMT not only mod- 
els spectral statistical fluctuation properties in chaotic 
quantum systems, but goes beyond and applies as well 
to chaotic classical systems. 

(ii) The distribution of length matrix elements P{A) it- 
self is universal to leading order. The difference be- 
tween the shape of P{A) and a Gaussian distribution 
gives a quantitative measure of how much the motion 
of the billiard ball in a chaotic billiard differs from a 
random walk. 

(iii) For the integrable rectangular billiard we find in 
correlation coefficient C 1 , in spectral rigidity A3 (L) 
and in the NNS distribution P{s) strong evidence for 
rigid behavior and strong correlation between neighbor 
eigenvalues. Such behavior can be understood from the 
observation that integrability introduces strong correla- 
tions in length matrix elements. This proliferates to the 
spectra. 

(iv) In contrast, for integrable quantum systems the 
NNS distribution generally shows Poissonian behavior 
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with correlation coefficient C = (there are a few ex- 
eptions). Thus from the point of view of level spacing 
distributions, a marked difference shows up between the 
classical and the quantum world. In our opinion, such 
difference is due to the group properties/symmetries 
coming from the loction of the boundary points of clas- 
sical trajectories. 
Open questions: 

(i) We have asked the question if universality holds 
beyond billiard systems, in particular for classically 
chaotic potential systems. We have studied an opti- 
cal billiard having some similarity to a potential sys- 
tem. First tentative results hint to the possibility that 
universality holds. Is this an indicator of possible diffu- 
sive/random walk behavior underlying such system? 

(ii) All the work proving universality in quantum chaos 
relies heavily on Gutzwiller's semi-classical trace for- 
mula and in particular is based on periodic orbits ll9[ fTQ[ 

Il4l 11511 . Likewise all the work establishing 
universality in classical chaos via correlation functions 
of periodic orbits OHIIllIlSHIIilli,!^ 
is based on periodic orbits. In contrast, our study of 
chaotic classical billiards, universality is observed based 



on non-periodic bouncing trajectories. How can finger- 
prints of chaos show up in periodic orbits as well as in 
non-periodic bouncing orbits? 
Future directions: 

(i) We plan to do numerical studies of the periodic 
Lorentz gas with finite horizon, i.e. in the regime where 
mathematically rigorous results have been established. 
We want to see if the character of length distribution 
changes in such regime. 

(ii) We hope that these findings may contribute to ob- 
tain a unified description of both, quantum and classical 
chaos, and help understanding why quantum chaos is 
typically weaker than classical chaos, e.g., via an effec- 
tive quantum action Js^, Ell] . 

(iii) The global statistical approach to classical chaos 
proposed here may help to give insight into the prob- 
lem of ergodicity breaking in Hamiltonian systems (e.g., 
dense packing of discs in the Lorentz gas model 1I22I1 ). 
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